#include "stdafx.h"
#include "OvTestFrSmart.h"
#include "ovbase/ovdebug.h"
#include "ovbase/ovfileio.h"
#include "ovbase/ovalgorithm.h"
#include "ovbase/ovproperty.h"
#include "ovutils/gridmapping.h"
#include "ovutils/ovfileconverter.h"
#include "ovutils/cpgresize.h"
#include "ovutils/cpgboundary.h"

namespace OVTEST_CASE {
    int OvTestFrSmart::RunQ1(int argc, char* argv[])
    {
        QString main_dir;
        QString gridf_name;
        QString gridc_name;
        QString propf_name;
        QString propc_name;
        if (argc < 7) {
            main_dir = QStringLiteral(R"(D:\PETROLEUM\RnD\FrSmart\Question1\)");
            gridf_name = "Grid_A.out";
            gridc_name = "Grid_B.out";
            propf_name = "A_porps.txt";
            propc_name = "B_porps.txt";
        }
        else {
            main_dir   = argv[2];
            gridf_name = argv[3];
            propf_name = argv[4];
            gridc_name = argv[5];
            propc_name = argv[6];
        }

        OVDEBUG = std::make_shared<OvDebug>();
        OVDEBUG->Setup("NULL", "OVTEST_CASE", NOV::DEBUG_STD, 5, 5, false);
        OvDebugStart("FrSmart_Q1");

        OvGrdecl gridf, gridc;
        OvPropList propsf, propsc;
        bool finish_read = false;
        if (!finish_read) {
            //OvEclIO::SetupGrdecl(FullFileName(main_dir, gridf_name), &gridf, nullptr);
            //OvEclIO::SetupGrdecl(FullFileName(main_dir, gridc_name), &gridc, nullptr);
            //OvEclIO::SetupGrdeclNew(FullFileName(main_dir, gridf_name), &gridf, nullptr);
            //OvEclIO::SetupGrdeclNew(FullFileName(main_dir, gridc_name), &gridc, nullptr);
            OvEclIO::ReadEclDataFile(FullFileName(main_dir, gridf_name), gridf, nullptr, QString(""), true);
            OvEclIO::ReadEclDataFile(FullFileName(main_dir, gridc_name), gridc, nullptr, QString(""), true);
            OvEclIO::ReadEclStylePropFile(FullFileName(main_dir, propf_name), &propsf, gridf.NI(), gridf.NJ(), gridf.NK());
            //OvEclIO::SaveGrdecl(&gridf, &propsf, FullFileName(main_dir, "A.GRDECL"));
        }
        else {
            gridf.Load(FullFileName(main_dir, "gf.GRDECL.B"));
        }

        OvPoints censf = gridf.GetCenterPoints();
        IMap im;
        GridMapping gm;
        gm.SetGridX(&censf, 0);
        gm.SetGridY(&gridc);
        gm.Map2(im);

        {
            const std::vector<int>& actnumf = gridf.Actnum();
            const std::vector<int>& actnumc = gridc.Actnum();
            const int sizec = gridc.Size();
            QStringList names = propsf.NamesSorted();
            for (auto name : names) {
                OvProp* prop = propsf.Property(name);
                if (prop->Type() != NOV::VARIANT_DOUBLE) continue;
                QVector<double> sum(sizec, 0.0);
                QVector<int> num(sizec, 0);
                std::vector<double> data = prop->DataVecD();
                for (int f = 0; f < data.size(); f++) {
                    int c = im[f];
                    if (actnumc[c] > 0 && actnumf[f] > 0) {
                        sum[c] += data[f];
                        num[c]++;
                    }
                }
                for (int c = 0; c < sizec; c++) {
                    if (num[c] > 0) {
                        sum[c] /= num[c];
                    }
                }
                propsc.AddMemoryCopy(name, sum.data(), sum.size(), sum.size() * sizeof(double), NOV::VARIANT_DOUBLE);
            }
        }

        //OvPropList propsc;
        //CpgCoarsening cpgc;
        //cpgc.SetGrid(&gridf, &propsf);
        //cpgc.SetCoarseScale(1, 1, 40);
        //cpgc.Coarsening(&gridc, &propsc);
        propsc.SaveEclipseKeywordsSingleFile(FullFileName(main_dir, propc_name));
        //gridc.Save(FullFileName(main_dir, "gridc.GRDECL.B"));
        //gridc.SaveAscii(FullFileName(main_dir, "gridc.GRDECL"));

        OvDebugFinish();
        OVDEBUG = nullptr;
        return 0;
    }
    int OvTestFrSmart::RunQ2(int argc, char* argv[])
    {
        QString main_dir;
        QString iname;
        QString ogname;
        QString opname;
        if (argc < 6) {
            main_dir = QStringLiteral(R"(D:\PETROLEUM\RnD\FrSmart\Question2\)");
            iname  = "Grid_and_data.txt";
            ogname = "Result_grid_coordinates.txt";
            opname = "Result_grid_S1.txt";
        }
        else {
            main_dir = argv[2];
            iname    = argv[3];
            ogname   = argv[4];
            opname   = argv[5];
        }

        OVDEBUG = std::make_shared<OvDebug>();
        OVDEBUG->Setup("NULL", "OVTEST_CASE", NOV::DEBUG_STD, 5, 5, false);
        OvDebugStart("FrSmart_Q2");

        OvPoints ptsi, ptso;
        std::vector<double> s1i, s1o;

        OvDebugStart("Read input data");
        Q2ReadGrid1(FullFileName(main_dir, iname), ptsi, s1i);
        OvDebugFinish();

        OvDebugStart("Read output data");
        Q2ReadGrid2(FullFileName(main_dir, ogname), ptso);
        OvDebugFinish();

        OvDebugStart("Interpolating S1");
        int mode = 2;
        if (mode == 1) {
            IMap im;
            GridMapping gm;
            gm.SetGridX(&ptsi, 0);
            gm.SetGridY(&ptso);
            gm.Map2(im);
            {
                const int sizec = ptso.size();
                s1o.assign(sizec, 0.0);
                std::vector<int> num(sizec, 0);
                for (int f = 0; f < s1i.size(); f++) {
                    int c = im[f];
                    s1o[c] += s1i[f];
                    num[c]++;
                }
                for (int c = 0; c < sizec; c++) {
                    if (num[c] > 0) {
                        s1o[c] /= num[c];
                    }
                }
            }
        }
        else {
            QVector<QMap<double, int>> ims;
            int num = 6;
            Pts2PtsMapping(ptsi, ptso, num, ims);
            {
                const int sizec = ptso.size();
                s1o.assign(sizec, 0.0);
                for (int c = 0; c < sizec; c++) {
                    QMap<double, int> im = ims[c];
                    double wtsum = 0.0;
                    for (auto it = im.begin(); it != im.end(); it++) {
                        double dist = it.key();
                        int f = it.value();
                        if (dist < 1e-3) {
                            s1o[c] = s1i[f];
                            break;
                        }
                        else {
                            double wt = 1.0 / dist / dist;
                            s1o[c] += s1i[f] * wt;
                            wtsum += wt;
                        }
                    }
                    if (wtsum > 1e-10) {
                        s1o[c] /= wtsum;
                    }
                }
            }
        }
        OvDebugFinish();

        OvDebugStart("Write output data");
        Q2WriteGrid2(FullFileName(main_dir, opname), ptso, s1o);
        OvDebugFinish();

        OvDebugFinish();
        OVDEBUG = nullptr;
        return 0;
    }

    int OvTestFrSmart::RunQ31(int argc, char* argv[])
    {
        QString main_dir;
        QString cpg_name;
        QString fem_name = "finiteGrid_Hex.vtk";
        if (argc < 4) {
            main_dir = QStringLiteral(R"(D:\PETROLEUM\RnD\FrSmart\Question3\)");
            cpg_name = "Grid_E.out";
        }
        else {
            main_dir = argv[2];
            cpg_name = argv[3];
        }

        OVDEBUG = std::make_shared<OvDebug>();
        OVDEBUG->Setup("NULL", "OVTEST_CASE", NOV::DEBUG_STD, 2, 2, false);
        OvDebugStart("FrSmart_Q3");

        OvDebugStart("Read input data");
        OvGrdecl cpg;
        OvEclIO::ReadEclDataFile(FullFileName(main_dir, cpg_name), cpg, nullptr, QString(""), true);
        OvDebugFinish();

        OvDebugStart("Analysis faults");
        OvPoints nodes;
        GetFaults1(cpg, nodes);
        OvDebugFinish();

        OvDebugStart("Writing VTK file");
        CpgSaveVtkFile(cpg, FullFileName(main_dir, fem_name), nodes);
        OvDebugFinish();

        OVDEBUG = nullptr;
        return 0;
    }


    int OvTestFrSmart::RunQ32(int argc, char* argv[])
    {
        QString main_dir;
        QString cpg_name;
        QString fem_name = "finiteGrid_Tet.vtk";
        if (argc < 4) {
            main_dir = QStringLiteral(R"(D:\PETROLEUM\RnD\FrSmart\Question3\)");
            cpg_name = "Grid_E.out";
        }
        else {
            main_dir = argv[2];
            cpg_name = argv[3];
        }

        OVDEBUG = std::make_shared<OvDebug>();
        OVDEBUG->Setup("NULL", "OVTEST_CASE", NOV::DEBUG_STD, 2, 2, false);
        OvDebugStart("FrSmart_Q3");

        OvDebugStart("Read input data");
        OvGrdecl cpg;
        OvEclIO::ReadEclDataFile(FullFileName(main_dir, cpg_name), cpg, nullptr, QString(""), true);
        OvDebugFinish();

        OvDebugStart("Analysis faults");
        OvPoints nodes;
        OvFaces faces;
        GetFaults2(cpg, nodes, faces);
        OvFracture3D frac;
        frac.SetData2(nodes, faces, nullptr, nullptr, nullptr, nullptr);
        frac.SetMarker();
        frac.Shrink(-0.01);
        OvDebugFinish();

        OvDebugStart("Generating Tetrahedral grids");
        QString meshing_command = "pqfeenna3000";
        //QString meshing_command = "pqfeenn";
        OvMeshPtr mesh = Q3MeshingTet(main_dir, cpg, frac, OvPropList(), meshing_command);
        OvDebugFinish();

        OvDebugStart("Writing VTK file");
        mesh->SaveVtkFile(FullFileName(main_dir, fem_name));
        OvDebugFinish();

        OVDEBUG = nullptr;
        return 0;
    }





























    void OvTestFrSmart::Q2ReadGrid1(QString fn, OvPoints& pts, std::vector<double>& s1)
    {
        pts.clear();
        s1.clear();
        //QRegExp delimiter(REGEXP_BLANK);
        QRegExp delimiter("[ \t]");
        //QChar delimiter('\t');
        QFile fdata(fn);
        if (!fdata.open(QIODevice::ReadOnly)) {
            OvDebugError("Cannot open " + fn, 1);
        }
        OvDebugInfo(QString("File Size: %1 mb").arg(fdata.size() / 1024.0 / 1024.0));
        QTextStream ts(&fdata);
        ts.readLine();
        while (!ts.atEnd()) {
            QString lll = ts.readLine();
            QString buff = lll.trimmed();
            if (buff.isEmpty()) continue;
            QStringList bufflist = buff.split(delimiter, QString::SkipEmptyParts);
            if (!bufflist.size() == 4) continue;
            double x = bufflist.at(0).toDouble();
            double y = bufflist.at(1).toDouble();
            double z = bufflist.at(2).toDouble();
            double s = bufflist.at(3).toDouble();
            pts.push_back(OvPoint(x, y, z));
            s1.push_back(s);
        }
        OvDebugInfo(QString("Point num = %1").arg(pts.size()));
        fdata.close();
    }

    void OvTestFrSmart::Q2ReadGrid2(QString fn, OvPoints& pts)
    {
        pts.clear();
        //QRegExp delimiter(REGEXP_BLANK);
        QRegExp delimiter("[ \t]");
        //QChar delimiter('\t');
        QFile fdata(fn);
        if (!fdata.open(QIODevice::ReadOnly)) {
            OvDebugError("Cannot open " + fn, 1);
        }
        OvDebugInfo(QString("File Size: %1 mb").arg(fdata.size() / 1024.0 / 1024.0));
        QTextStream ts(&fdata);
        ts.readLine();
        while (!ts.atEnd()) {
            QString lll = ts.readLine();
            QString buff = lll.trimmed();
            if (buff.isEmpty()) continue;
            QStringList bufflist = buff.split(delimiter, QString::SkipEmptyParts);
            if (!bufflist.size() == 4) continue;
            double x = bufflist.at(1).toDouble();
            double y = bufflist.at(2).toDouble();
            double z = bufflist.at(3).toDouble();
            pts.push_back(OvPoint(x, y, z));
        }
        fdata.close();
        OvDebugInfo(QString("Point num = %1").arg(pts.size()));
    }

    void OvTestFrSmart::Q2WriteGrid2(QString fn, const OvPoints& pts, const std::vector<double>& s1)
    {
        const QChar delimiter('\t');
        QFile fdata(fn);
        if (!fdata.open(QIODevice::WriteOnly)) {
            OvDebugError("Cannot open " + fn, 1);
        }
        QTextStream ts(&fdata);
        ts << "x	y	z	S1\n";
        for (int i = 0; i < pts.size(); i++) {
            ts << pts.at(i).X() << delimiter
                << pts.at(i).Y() << delimiter
                << pts.at(i).Z() << delimiter
                << s1.at(i) << '\n';
        }
        fdata.close();
        OvDebugInfo(QString("Point num = %1").arg(pts.size()));
    }

    void OvTestFrSmart::Pts2PtsMapping(const OvPoints& fcens, const OvPoints& ccens, int num, QVector<QMap<double, int>>& ims)
    {
        const int fsize = fcens.size();
        const int csize = ccens.size();
        std::vector<double> flens(fsize);
        for (int i = 0; i < fsize; i++) {
            flens[i] = fcens[i].Length();
        }
        std::vector<int> orders;
        NOV::sort_indexes_2(flens, orders);
        OvPoints fcens_orders(fsize);
        for (int i = 0; i < fsize; i++) {
            fcens_orders[i] = fcens[orders[i]];
        }
        ims.resize(csize);
        const int max_thread_num = std::thread::hardware_concurrency();
        const int use_thread_num = std::min(int(max_thread_num * 0.75), csize / 100);
        OvDebugInfo(QString("point_num = %1, use_thread_num=%2").arg(csize).arg(use_thread_num));
        int percent_num = 5;
#pragma omp parallel for num_threads(use_thread_num)
        for (int i = 0; i < csize; i++) {
            OvPoint pt = ccens[i];
            QMap<double, int> im;
            NOV::FindNearestPoints(fsize, fcens_orders.data(), flens.data(), pt, num, im);
            for (QMap<double, int>::iterator it = im.begin(); it != im.end(); it++) {
                it.value() = orders.at(it.value());
            }
            ims[i] = im;
        }
    }





    OvMeshPtr OvTestFrSmart::Q3MeshingTet(QString main_dir, OvGrdecl cpg, OvFracture3D frac, OvPropList prop_list, QString meshing_command)
    {
        //const QString mesh_dir = Dir(main_dir, "MESH");
        const QString model_name = "tet";

        QVector<OvFracture3D> fracs;
        QVector<OvWfMeshing::FractureType> frac_types;
        fracs.append(frac);
        frac_types.append(OvWfMeshing::FRAC_DISCRETE);
        //for (int i = 0; i < frac_sets_.size(); i++) {
        //    FractureSet& fs = frac_sets_.at(i);
        //    fracs.append(*fs.frac);
        //    frac_types.append(fs.type);
        //}

        OvBoundary bd;
        CpgBoundary cb;
        cb.SetGrid(&cpg);
        cb.GetBoundary3D(bd, true, false, 1000);

        //QString mat_ts_name = FullFileName(mesh_dir, "matrix.ts");
        //bd.SaveGocadTsFile(mat_ts_name);

        bool real_edfm = false;
        OvWfMeshing meshing(&fracs, &frac_types, real_edfm, &bd, meshing_command, main_dir, model_name, &cpg, false);
        meshing.Start();
        OvMeshPtr mesh = meshing.Mesh();
        //OvMesher::Write3DMesh2TecplotGeneral(mesh.get(), nullptr, false, FullFileName(mesh_dir, model_name + ".dat"));
        //mesh->SaveAscii(FullFileName(mesh_dir, model_name + MESH_SUFFIX_A));
        //mesh->Save(FullFileName(mesh_dir, model_name + MESH_SUFFIX_B));
        return mesh;
    }

    void OvTestFrSmart::GetFaults1(OvGrdecl cpg, OvPoints& nodes)
    {
        cpg.SmoothGridsByActnum();
        cpg.LoadNodes();
        const std::vector<int>& actnum = cpg.Actnum();
        const int ni = cpg.NI();
        const int nj = cpg.NJ();
        std::vector<double>& zcorn = cpg.Zcorn();
        nodes.clear();
        for (int n = 0; n < cpg.Size(); n++) {
            size_t i, j, k;
            cpg.GetIJK(n, i, j, k);
            OvPoints hex(8);
            cpg.CalcNodes(n, hex.data());
            std::vector<int> top(4);
            cpg.NodeIndexOnBoundaryFace(NOV::POS_TOP, top.data());
            OvPolygon plg1, plg2;
            plg1.append(hex[top.at(0)]);
            plg1.append(hex[top.at(1)]);
            plg1.append(hex[top.at(3)]);
            plg1.append(hex[top.at(2)]);
            plg2.append(hex[top.at(0) + 4]);
            plg2.append(hex[top.at(1) + 4]);
            plg2.append(hex[top.at(3) + 4]);
            plg2.append(hex[top.at(2) + 4]);
            OvLine e;
            int idegen = -1;
            std::vector<int> ei;
            bool flag = plg1.FindDegeneratedVertex(ei, idegen, e);
            if (flag) {
                nodes.push_back(plg1.Node(idegen));
                nodes.push_back(plg2.Node(idegen));
            }
        }
    }

    void OvTestFrSmart::GetFaults2(OvGrdecl cpg, OvPoints& nodes, OvFaces& faces)
    {
        cpg.SmoothGridsByActnum();
        cpg.LoadNodes();

        const std::vector<int>& actnum = cpg.Actnum();
        const int ni = cpg.NI();
        const int nj = cpg.NJ();
        nodes.clear();
        faces.clear();
        const int size = cpg.Size();
        OvPosilist pos = cpg.GetPosiList();
        for (int idx = 0; idx < size; idx++) {
            if (!actnum.at(idx) || (!pos.at(idx).Contains(NOV::POS_TOP) && !pos.at(idx).Contains(NOV::POS_BOTTOM)))
                continue;
            OvPoint hex[8];
            cpg.CalcNodes(idx, hex);
            for (int i = 0; i < 6; i++) {
                NOV::EPosition p = cpg.BoundaryPosition(i);
                if (p != NOV::POS_NORTH && p != NOV::POS_SOUTH) continue;
                if (!pos.at(idx).Contains(p)) continue;

                            
                if (pos.at(idx).Contains(NOV::POS_TOP)) {
                    std::vector<int> top(4);
                    cpg.NodeIndexOnBoundaryFace(NOV::POS_TOP, top.data());
                    OvPolygon plg;  
                    plg.append(hex[top.at(0)]);
                    plg.append(hex[top.at(1)]);
                    plg.append(hex[top.at(3)]);
                    plg.append(hex[top.at(2)]);
                    OvLine e;
                    int idegen = -1;
                    std::vector<int> ei;
                    bool flag = plg.FindDegeneratedVertex(ei, idegen, e);
                    if (!flag) continue;
                    if (p == NOV::POS_NORTH && (idegen == 2 || idegen == 3)) {
                        OvPoint pt = plg.Node(idegen);
                        nodes.push_back(pt);
                    }
                    if (p == NOV::POS_SOUTH && (idegen == 0 || idegen == 1)) {
                        OvPoint pt = plg.Node(idegen);
                        nodes.push_back(pt);
                    }
                }

                if (pos.at(idx).Contains(NOV::POS_BOTTOM)) {
                    std::vector<int> top(4);
                    cpg.NodeIndexOnBoundaryFace(NOV::POS_BOTTOM, top.data());
                    OvPolygon plg;
                    plg.append(hex[top.at(0)]);
                    plg.append(hex[top.at(1)]);
                    plg.append(hex[top.at(3)]);
                    plg.append(hex[top.at(2)]);
                    OvLine e;
                    int idegen = -1;
                    std::vector<int> ei;
                    bool flag = plg.FindDegeneratedVertex(ei, idegen, e);
                    if (!flag) continue;
                    if (p == NOV::POS_NORTH && (idegen == 2 || idegen == 3)) {
                        OvPoint pt = plg.Node(idegen);
                        nodes.push_back(pt);
                    }
                    if (p == NOV::POS_SOUTH && (idegen == 0 || idegen == 1)) {
                        OvPoint pt = plg.Node(idegen);
                        nodes.push_back(pt);
                    }
                }
            }
        }

        if (nodes.size() == 4) {
            OvPolygon plg;
            plg.SetNodes(nodes.size(), nodes.data());
            plg.SortAnitclockwise();
            for (int i = 0; i < nodes.size(); i++) nodes[i] = plg.Node(i);
            faces.push_back(OvMeshObject(0, 1, 2, 3, 1));
        }
    }


    void OvTestFrSmart::CpgSaveVtkFile(OvGrdecl cpg, QString file_name, OvPoints& fault_nodes)
    {
        cpg.LoadNodes();
        OvPoints nodes;
        OvFaces faces;
        for (int n = 0; n < cpg.Size(); n++) {
            if (!cpg.CellActnum(n)) continue;
            OvPoints hex(8);
            cpg.CalcNodes(n, hex.data());
            int ibase = nodes.size();
            nodes.push_back(hex.at(0));
            nodes.push_back(hex.at(1));
            nodes.push_back(hex.at(3));
            nodes.push_back(hex.at(2));
            nodes.push_back(hex.at(4));
            nodes.push_back(hex.at(5));
            nodes.push_back(hex.at(7));
            nodes.push_back(hex.at(6));
            std::vector<int> l(8);
            for (int i = 0; i < l.size(); i++) l[i] = ibase + i;
            faces.push_back(OvMeshObject(l.size(), l.data(), 0));
        }
        QVector<int> point_map;
        NOV::RemoveDuplicatePointsInFaces(nodes, faces, point_map);
        
        IMap im;
        GridMapping gm;
        gm.SetGridX(&fault_nodes, 0);
        gm.SetGridY(&nodes);
        gm.Map2(im);

        const int node_num = nodes.size();
        const int cell_num = faces.size();
        QFile f(file_name);
        if (!f.open(QFile::WriteOnly))
            OvDebugError(QString("Cannot Write File: %1").arg(file_name), 1);
        QTextStream ts(&f);
        ts << QString("# vtk DataFile Version 2.0\n");
        ts << QString("Unstructured Grid\n");
        ts << QString("ASCII\n");
        ts << QString("DATASET UNSTRUCTURED_GRID\n");
        ts << QString("POINTS %1 float\n").arg(node_num);
        for (int i = 0; i < node_num; i++) {
            OvPoint pt = nodes.at(i);
            ts << QString("%1 %2 %3\n").arg(pt.X(), 0, 'f', 3).arg(pt.Y(), 0, 'f', 3).arg(pt.Z(), 0, 'f', 3);
        }
        ts << "\n";
        ts << QString("CELLS %1 %2\n").arg(cell_num).arg(cell_num * 9);
        for (int i = 0; i < cell_num; i++) {
            OvMeshObject c = faces.at(i);
            ts << c.size();
            for (int j = 0; j < c.size(); j++) ts << " " << c.at(j);
            ts << "\n";
        }
        ts << "\n";
        ts << QString("CELL_TYPES %1\n").arg(cell_num);
        for (int i = 0; i < cell_num; i++) {
            ts << QString("12\n");
        }
        ts << "\n";

        std::vector<int> fault_flag(node_num, 0);
        for (int i = 0; i < im.size(); i++) {
            int element = im.at(i);
            fault_flag[element] = 1;
        }
        ts << QString("POINT_DATA %1\n").arg(node_num);
        ts << QString("SCALARS FAULT float 1\n");
        ts << QString("LOOKUP_TABLE default\n");
        for (int i = 0; i < node_num; i++) {
            ts << QString("%1\n").arg(fault_flag.at(i));
        }
        f.close();
    }
}